last updated: 2023-08-22

Getting Things Setup

As usual, make sure we have the right packages for this exercise

if (!require("pacman")) install.packages("pacman"); library(pacman)
## Loading required package: pacman
# let's load all of the files we were using and want to have again today
p_load("tidyverse", "knitr", "readr",
       "pander", "BiocManager", 
       "dplyr", "stringr", 
       "statmod", # required dependency, need to load manually on some macOS versions.
       "purrr", # for working with lists (beautify column names)
       "reactable") # for pretty tables.

# We also need these Bioconductor packages today.
p_load("edgeR", "AnnotationDbi", "org.Sc.sgd.db")
#NOTE: edgeR loads limma as a dependency

Description

This will be our last differential expression analysis workflow, converting gene counts across samples into meaningful information about genes that appear to be significantly differentially expressed between samples

Learning outcomes

At the end of this exercise, you should be able to:

library(limma)
library(org.Sc.sgd.db)
# for ease of use, set max number of digits after decimal
options(digits=3)

Loading in the featureCounts object

We saved this file in the last exercise (Read_Counting.Rmd) from the RSubread package. Now we can load that object back in and assign it to the variable fc. Be sure to change the file path if you have saved it in a different location.

path_fc_object <- path.expand("~/Desktop/Genomic_Data_Analysis/Data/yeast_fc_output.Rds")
fc <- readRDS(file = path_fc_object)

If you don’t have that file for any reason, the below code chunk will load a copy of it from Github.

if( !exists("fc") )
{
  fc <- read_rds('https://github.com/clstacy/GenomicDataAnalysis_Fa23/raw/main/data/ethanol_stress/yeast_fc_output.Rds')
}

To find the order of files we need, we can get just the part of the column name before the first “.” symbol with this command:

str_split_fixed(fc$counts %>% colnames(), "\\.", n = 2)[, 1]
##  [1] "YPS606_MSN24_ETOH_REP1_R1" "YPS606_MSN24_ETOH_REP2_R1"
##  [3] "YPS606_MSN24_ETOH_REP3_R1" "YPS606_MSN24_ETOH_REP4_R1"
##  [5] "YPS606_MSN24_MOCK_REP1_R1" "YPS606_MSN24_MOCK_REP2_R1"
##  [7] "YPS606_MSN24_MOCK_REP3_R1" "YPS606_MSN24_MOCK_REP4_R1"
##  [9] "YPS606_WT_ETOH_REP1_R1"    "YPS606_WT_ETOH_REP2_R1"   
## [11] "YPS606_WT_ETOH_REP3_R1"    "YPS606_WT_ETOH_REP4_R1"   
## [13] "YPS606_WT_MOCK_REP1_R1"    "YPS606_WT_MOCK_REP2_R1"   
## [15] "YPS606_WT_MOCK_REP3_R1"    "YPS606_WT_MOCK_REP4_R1"
sample_metadata <- tribble(
  ~Sample,                      ~Genotype,    ~Condition,
  "YPS606_MSN24_ETOH_REP1_R1",   "msn24dd",   "EtOH",
  "YPS606_MSN24_ETOH_REP2_R1",   "msn24dd",   "EtOH",
  "YPS606_MSN24_ETOH_REP3_R1",   "msn24dd",   "EtOH",
  "YPS606_MSN24_ETOH_REP4_R1",   "msn24dd",   "EtOH",
  "YPS606_MSN24_MOCK_REP1_R1",   "msn24dd",   "unstressed",
  "YPS606_MSN24_MOCK_REP2_R1",   "msn24dd",   "unstressed",
  "YPS606_MSN24_MOCK_REP3_R1",   "msn24dd",   "unstressed",
  "YPS606_MSN24_MOCK_REP4_R1",   "msn24dd",   "unstressed",
  "YPS606_WT_ETOH_REP1_R1",      "WT",        "EtOH",
  "YPS606_WT_ETOH_REP2_R1",      "WT",        "EtOH",
  "YPS606_WT_ETOH_REP3_R1",      "WT",        "EtOH",
  "YPS606_WT_ETOH_REP4_R1",      "WT",        "EtOH",
  "YPS606_WT_MOCK_REP1_R1",      "WT",        "unstressed",
  "YPS606_WT_MOCK_REP2_R1",      "WT",        "unstressed",
  "YPS606_WT_MOCK_REP3_R1",      "WT",        "unstressed",
  "YPS606_WT_MOCK_REP4_R1",      "WT",        "unstressed") %>%
  # Create a new column that combines the Genotype and Condition value
  mutate(Group = factor(
    paste(Genotype, Condition, sep = "."),
    levels = c(
      "WT.unstressed","WT.EtOH",
      "msn24dd.unstressed", "msn24dd.EtOH"
    )
  )) %>%
  # make Condition and Genotype a factor (with baseline as first level) for edgeR
  mutate(
    Genotype = factor(Genotype,
                      levels = c("WT", "msn24dd")),
    Condition = factor(Condition,
                       levels = c("unstressed", "EtOH"))
  )

Now, let’s create a design matrix with this information

group <- sample_metadata$Group
design <- model.matrix(~ 0 + group)

# beautify column names
colnames(design) <- levels(group)
design
##    WT.unstressed WT.EtOH msn24dd.unstressed msn24dd.EtOH
## 1              0       0                  0            1
## 2              0       0                  0            1
## 3              0       0                  0            1
## 4              0       0                  0            1
## 5              0       0                  1            0
## 6              0       0                  1            0
## 7              0       0                  1            0
## 8              0       0                  1            0
## 9              0       1                  0            0
## 10             0       1                  0            0
## 11             0       1                  0            0
## 12             0       1                  0            0
## 13             1       0                  0            0
## 14             1       0                  0            0
## 15             1       0                  0            0
## 16             1       0                  0            0
## attr(,"assign")
## [1] 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"

Count loading and Annotation

The count matrix is used to construct a DGEList class object. This is the main data class in the edgeR package. The DGEList object is used to store all the information required to fit a generalized linear model to the data, including library sizes and dispersion estimates as well as counts for each gene.

y <- DGEList(fc$counts, group=group)
colnames(y) <- sample_metadata$Sample
y$samples
##                                        group lib.size norm.factors
## YPS606_MSN24_ETOH_REP1_R1       msn24dd.EtOH   195105            1
## YPS606_MSN24_ETOH_REP2_R1       msn24dd.EtOH   180513            1
## YPS606_MSN24_ETOH_REP3_R1       msn24dd.EtOH   165696            1
## YPS606_MSN24_ETOH_REP4_R1       msn24dd.EtOH   172157            1
## YPS606_MSN24_MOCK_REP1_R1 msn24dd.unstressed   137673            1
## YPS606_MSN24_MOCK_REP2_R1 msn24dd.unstressed   141652            1
## YPS606_MSN24_MOCK_REP3_R1 msn24dd.unstressed   171295            1
## YPS606_MSN24_MOCK_REP4_R1 msn24dd.unstressed   172884            1
## YPS606_WT_ETOH_REP1_R1               WT.EtOH   154719            1
## YPS606_WT_ETOH_REP2_R1               WT.EtOH   169239            1
## YPS606_WT_ETOH_REP3_R1               WT.EtOH   179811            1
## YPS606_WT_ETOH_REP4_R1               WT.EtOH   157740            1
## YPS606_WT_MOCK_REP1_R1         WT.unstressed   186150            1
## YPS606_WT_MOCK_REP2_R1         WT.unstressed   154835            1
## YPS606_WT_MOCK_REP3_R1         WT.unstressed   184436            1
## YPS606_WT_MOCK_REP4_R1         WT.unstressed   172568            1

Human-readable gene symbols can also be added to complement the gene ID for each gene, using the annotation in the org.Sc.sgd.db package.

y$genes <- AnnotationDbi::select(org.Sc.sgd.db,keys=rownames(y),columns="GENENAME")
## 'select()' returned 1:1 mapping between keys and columns
head(y$genes)
##       ORF        SGD GENENAME
## 1 YDL246C S000002405     SOR2
## 2 YDL243C S000002402     AAD4
## 3 YDR387C S000002795    CIN10
## 4 YDL094C S000002252     <NA>
## 5 YDR438W S000002846    THI74
## 6 YDR523C S000002931     SPS1

Filtering to remove low counts

Genes with very low counts across all libraries provide little evidence for differential ex- pression. In addition, the pronounced discreteness of these counts interferes with some of the statistical approximations that are used later in the pipeline. These genes should be filtered out prior to further analysis. Here, we will retain a gene only if it is expressed at a count-per-million (CPM) above 60 in at least four samples.

keep <- rowSums(cpm(y) > 60) >= 4
y <- y[keep,]
summary(keep)
##    Mode   FALSE    TRUE 
## logical    4518    2609

Where did those cutoff numbers come from?

As a general rule, we don’t want to exclude a gene that is expressed in only one group, so a cutoff number equal to the number of replicates can be a good starting point. For counts, a good threshold can be chosen by identifying the CPM that corresponds to a count of 10, which in this case would be about 60 (due to our fastq files being subsets of the full reads):

cpm(10, mean(y$samples$lib.size))
##      [,1]
## [1,] 59.3

Smaller CPM thresholds are usually appropriate for larger libraries.

Normalization for composition bias

TMM normalization is performed to eliminate composition biases between libraries. This generates a set of normalization factors, where the product of these factors and the library sizes defines the effective library size. The calcNormFactors function returns the DGEList argument with only the norm.factors changed.

y <- calcNormFactors(y)
y$samples
##                                        group lib.size norm.factors
## YPS606_MSN24_ETOH_REP1_R1       msn24dd.EtOH   195105        1.156
## YPS606_MSN24_ETOH_REP2_R1       msn24dd.EtOH   180513        1.094
## YPS606_MSN24_ETOH_REP3_R1       msn24dd.EtOH   165696        1.074
## YPS606_MSN24_ETOH_REP4_R1       msn24dd.EtOH   172157        0.997
## YPS606_MSN24_MOCK_REP1_R1 msn24dd.unstressed   137673        1.038
## YPS606_MSN24_MOCK_REP2_R1 msn24dd.unstressed   141652        1.046
## YPS606_MSN24_MOCK_REP3_R1 msn24dd.unstressed   171295        0.999
## YPS606_MSN24_MOCK_REP4_R1 msn24dd.unstressed   172884        0.996
## YPS606_WT_ETOH_REP1_R1               WT.EtOH   154719        0.865
## YPS606_WT_ETOH_REP2_R1               WT.EtOH   169239        0.908
## YPS606_WT_ETOH_REP3_R1               WT.EtOH   179811        0.945
## YPS606_WT_ETOH_REP4_R1               WT.EtOH   157740        0.928
## YPS606_WT_MOCK_REP1_R1         WT.unstressed   186150        1.004
## YPS606_WT_MOCK_REP2_R1         WT.unstressed   154835        1.065
## YPS606_WT_MOCK_REP3_R1         WT.unstressed   184436        0.942
## YPS606_WT_MOCK_REP4_R1         WT.unstressed   172568        0.984

The normalization factors multiply to unity across all libraries. A normalization factor below unity indicates that the library size will be scaled down, as there is more suppression (i.e., composition bias) in that library relative to the other libraries. This is also equivalent to scaling the counts upwards in that sample. Conversely, a factor above unity scales up the library size and is equivalent to downscaling the counts. The performance of the TMM normalization procedure can be examined using mean- difference (MD) plots. This visualizes the library size-adjusted log-fold change between two libraries (the difference) against the average log-expression across those libraries (the mean). The below command plots an MD plot, comparing sample 1 against an artificial library constructed from the average of all other samples.

for (sample in 1:nrow(y$samples)) {
  plotMD(cpm(y, log=TRUE), column=sample)
  abline(h=0, col="red", lty=2, lwd=2)
}

Exploring differences between libraries

The data can be explored by generating multi-dimensional scaling (MDS) plots. This visualizes the differences between the expression profiles of different samples in two dimensions. The next plot shows the MDS plot for the yeast heatshock data.

points <- c(1,1,2,2)
colors <- rep(c("black", "red"),8)
plotMDS(y, col=colors[group], pch=points[group])
legend("topright", legend=levels(group),
     pch=points, col=colors, ncol=2)

Estimate Dispersion

This is the first step in a limma analysis that differs from the edgeR workflow.

y <- voom(y, design, plot = T)

# compare this to the edgeR function estimateDisp, which uses a NB distribution.
# y <- estimateDisp(y, design, robust=TRUE)
# plotBCV(y)

What is voom doing?

Limma uses the lmFit function. This returns a MArrayLM object containing the weighted least squares estimates for each gene.

fit <- lmFit(y, design)
head(coef(fit))
##         WT.unstressed WT.EtOH msn24dd.unstressed msn24dd.EtOH
## YDR492W          6.77    3.28               6.73         2.76
## YDR508C          9.13    8.59               9.02         9.01
## YDR186C          5.63    6.53               5.78         6.38
## YDR150W          7.77    6.74               7.87         6.57
## YDL182W          8.96    6.63               8.70         7.31
## YDR232W          8.03    7.56               7.68         7.07
# edgeR equivalent
# fit <- glmQLFit(y, design, robust=TRUE)
# head(fit$coefficients)
# plotQLDisp(fit)

Comparisons between groups (log fold-changes) are obtained as contrasts of these fitted linear models:

Testing for differential expression

The final step is to actually test for significant differential expression in each gene, using the QL F-test. The contrast of interest can be specified using the makeContrasts function in limma, the same one that is used by edgeR.

# generate contrasts we are interested in learning about
my.contrasts <- makeContrasts(EtOHvsMOCK.WT = WT.EtOH - WT.unstressed, 
                     EtOHvsMOCK.MSN24dd = msn24dd.EtOH - msn24dd.unstressed,
                     EtOH.MSN24ddvsWT = msn24dd.EtOH - WT.EtOH,
                     MOCK.MSN24ddvsWT = msn24dd.unstressed - WT.unstressed,
                     EtOHvsWT.MSN24ddvsWT = (msn24dd.EtOH-msn24dd.unstressed)-(WT.EtOH-WT.unstressed),
                     levels=design)

# fit the linear model to these contrasts
res_all <- contrasts.fit(fit, my.contrasts)

# This looks at all of our contrasts in my.contrasts
res_all <- eBayes(res_all)

# eBayes is the alternative to glmQLFTest in edgeR
# This contrast looks at the difference in the stress responses between mutant and WT
# res <- glmQLFTest(fit, contrast = my.contrasts)
top.table <- topTable(res_all, sort.by = "F", n = Inf)
head(top.table, 20)
##             ORF        SGD GENENAME EtOHvsMOCK.WT EtOHvsMOCK.MSN24dd
## YJL034W YJL034W S000003571     KAR2          3.72              3.299
## YJL052W YJL052W S000003588     TDH1          3.46              2.674
## YGR254W YGR254W S000003486     ENO1          6.11              5.490
## YKL035W YKL035W S000001518     UGP1          4.63              0.616
## YIL053W YIL053W S000001315     GPP1          2.79              3.024
## YJR009C YJR009C S000003769     TDH2          2.39              2.136
## YBR126C YBR126C S000000330     TPS1          5.46              2.547
## YLR249W YLR249W S000004239     YEF3         -2.34             -1.930
## YIL169C YIL169C S000001431     CSS1         -2.20             -2.264
## YCR012W YCR012W S000000605     PGK1          2.17              1.773
## YCL040W YCL040W S000000545     GLK1          8.41              7.900
## YEL071W YEL071W S000000797     DLD3          2.49              3.317
## YGL008C YGL008C S000002976     PMA1         -2.39             -1.989
## YAL005C YAL005C S000000004     SSA1          2.90              1.519
## YCL043C YCL043C S000000548     PDI1          2.22              1.808
## YGL055W YGL055W S000003023     OLE1         -4.28             -3.666
## YMR217W YMR217W S000004830     GUA1         -3.38             -3.296
## YNL209W YNL209W S000005153     SSB2         -2.00             -1.681
## YPL265W YPL265W S000006186     DIP5          2.85              3.155
## YMR011W YMR011W S000004613     HXT2         -4.42             -4.627
##         EtOH.MSN24ddvsWT MOCK.MSN24ddvsWT EtOHvsWT.MSN24ddvsWT AveExpr    F
## YJL034W          -0.2220          0.19435              -0.4164   11.50 1000
## YJL052W          -0.8640         -0.07931              -0.7847   12.88  940
## YGR254W          -0.7536         -0.13617              -0.6174   11.35  759
## YKL035W          -4.0352         -0.02539              -4.0098    9.38  660
## YIL053W          -0.0631         -0.29550               0.2324   10.89  539
## YJR009C          -0.3295         -0.07751              -0.2520   13.74  528
## YBR126C          -3.4203         -0.50803              -2.9123    8.01  429
## YLR249W           0.3379         -0.06984               0.4077   12.43  415
## YIL169C           0.0432          0.10341              -0.0602   12.92  396
## YCR012W          -0.4585         -0.06682              -0.3917   13.10  353
## YCL040W          -2.0644         -1.55607              -0.5083    7.90  351
## YEL071W           0.6991         -0.12656               0.8257    9.40  338
## YGL008C           0.4125          0.00955               0.4029   12.28  335
## YAL005C          -0.9969          0.38262              -1.3795   11.42  329
## YCL043C          -0.4555         -0.04622              -0.4093   11.07  310
## YGL055W           0.5466         -0.06327               0.6099    8.79  303
## YMR217W           0.0289         -0.05075               0.0797    9.31  298
## YNL209W           0.2188         -0.09836               0.3171   12.05  297
## YPL265W           0.1944         -0.10685               0.3013    9.12  277
## YMR011W          -0.4367         -0.22653              -0.2102    8.73  262
##          P.Value adj.P.Val
## YJL034W 5.07e-50  1.32e-46
## YJL052W 2.98e-49  3.89e-46
## YGR254W 1.33e-46  1.15e-43
## YKL035W 6.87e-45  4.48e-42
## YIL053W 2.06e-42  1.08e-39
## YJR009C 3.65e-42  1.59e-39
## YBR126C 1.25e-39  4.67e-37
## YLR249W 3.08e-39  1.00e-36
## YIL169C 1.11e-38  3.20e-36
## YCR012W 2.78e-37  7.25e-35
## YCL040W 3.19e-37  7.56e-35
## YEL071W 9.21e-37  2.00e-34
## YGL008C 1.18e-36  2.37e-34
## YAL005C 1.81e-36  3.38e-34
## YCL043C 9.46e-36  1.65e-33
## YGL055W 1.84e-35  3.00e-33
## YMR217W 2.95e-35  4.29e-33
## YNL209W 2.96e-35  4.29e-33
## YPL265W 2.14e-34  2.94e-32
## YMR011W 9.34e-34  1.22e-31
top.table %>% 
  tibble() %>% 
  arrange(adj.P.Val) %>%
  mutate(across(where(is.numeric), signif, 3)) %>%
  reactable()
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(where(is.numeric), signif, 3)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))
# edgeR equivalent below:
# let's take a quick look at the results
# topTags(res, n=10) 
# 
# # generate a beautiful table for the pdf/html file.
# topTags(res, n=Inf) %>% data.frame() %>% 
#   arrange(FDR) %>%
#   mutate(logFC=round(logFC,2)) %>%
#   mutate(across(where(is.numeric), signif, 3)) %>%
#   reactable()
# Let's see how many genes in total are significantly different in any contrast
length(which(top.table$adj.P.Val < 0.05))
## [1] 1895
# let's summarize this and break it down by contrast.
res_all %>%
  decideTests(p.value = 0.05, lfc = 0) %>%
  summary()
##        EtOHvsMOCK.WT EtOHvsMOCK.MSN24dd EtOH.MSN24ddvsWT MOCK.MSN24ddvsWT
## Down             873                850              377                2
## NotSig           866               1016             2004             2606
## Up               870                743              228                1
##        EtOHvsWT.MSN24ddvsWT
## Down                    164
## NotSig                 2328
## Up                      117
# Bonus: limma allows us to create a venn diagram of these contrasts 
# up & downregulated genes
res_all %>%
  decideTests(p.value = 0.05, lfc = 1) %>% 
  vennDiagram(include=c("up", "down"),
              lwd=0.75,
              mar=rep(2,4), # increase margine size
              counts.col= c("red", "blue"),
              show.include=TRUE)

Looking at a specific contrast

It is interesting to see all of the contrasts simultaneously, but often we may want to look at just a single contrast (and get the corresponding probabilities). Here is how we do that:

# fit the linear model to these contrasts
res <- contrasts.fit(fit, my.contrasts[,"EtOHvsWT.MSN24ddvsWT"])

# This contrast looks at the difference in the stress responses between mutant and WT
res <- eBayes(res)
# Note that there is no longer an "F" column, because we only look at one contrast.
top.table <- topTable(res, sort.by = "P", n = Inf)
head(top.table, 20)
##             ORF        SGD GENENAME  logFC AveExpr      t  P.Value adj.P.Val
## YKL035W YKL035W S000001518     UGP1 -4.010    9.38 -19.98 9.76e-28  2.55e-24
## YPR149W YPR149W S000006353   NCE102 -3.857    7.58 -10.73 2.01e-15  2.63e-12
## YPL004C YPL004C S000005925     LSP1 -2.665    8.35 -10.59 3.36e-15  2.92e-12
## YDR343C YDR343C S000002751     HXT6 -5.274    7.30 -10.00 3.03e-14  1.87e-11
## YML100W YML100W S000004566     TSL1 -7.795    6.02  -9.95 3.59e-14  1.87e-11
## YDR077W YDR077W S000002484     SED1 -1.474   10.79  -9.79 6.42e-14  2.79e-11
## YMR105C YMR105C S000004711     PGM2 -6.545    6.22  -9.22 5.63e-13  1.82e-10
## YAL005C YAL005C S000000004     SSA1 -1.380   11.42  -9.19 6.21e-13  1.82e-10
## YKL150W YKL150W S000001633     MCR1 -3.068    7.77  -9.19 6.29e-13  1.82e-10
## YGR086C YGR086C S000003318     PIL1 -1.671    8.89  -8.78 2.94e-12  7.67e-10
## YBR126C YBR126C S000000330     TPS1 -2.912    8.01  -8.56 6.79e-12  1.61e-09
## YOL155C YOL155C S000005515     HPF1 -1.526   10.24  -8.52 7.95e-12  1.73e-09
## YJR045C YJR045C S000003806     SSC1 -1.243   10.38  -8.30 1.88e-11  3.77e-09
## YLL024C YLL024C S000003947     SSA2 -0.959   12.91  -7.81 1.22e-10  2.28e-08
## YDR133C YDR133C S000002540     <NA> -1.602    8.76  -7.72 1.77e-10  3.08e-08
## YFR053C YFR053C S000001949     HXK1 -7.158    4.54  -7.58 3.05e-10  4.98e-08
## YFL014W YFL014W S000001880    HSP12 -7.320    3.77  -7.41 5.94e-10  9.12e-08
## YMR196W YMR196W S000004809     <NA> -5.471    5.75  -7.37 6.84e-10  9.57e-08
## YJL078C YJL078C S000003614     PRY3 -1.674    8.90  -7.37 6.97e-10  9.57e-08
## YMR261C YMR261C S000004874     TPS3 -2.370    7.45  -7.35 7.35e-10  9.59e-08
##            B
## YKL035W 52.6
## YPR149W 24.6
## YPL004C 24.3
## YDR343C 21.6
## YML100W 19.5
## YDR077W 21.2
## YMR105C 17.7
## YAL005C 18.9
## YKL150W 19.2
## YGR086C 17.6
## YBR126C 16.8
## YOL155C 16.5
## YJR045C 15.6
## YLL024C 13.5
## YDR133C 13.6
## YFR053C 11.1
## YFL014W 10.5
## YMR196W 11.4
## YJL078C 12.2
## YMR261C 12.3
top.table %>% 
  tibble() %>% 
  arrange(adj.P.Val) %>%
  mutate(across(where(is.numeric), signif, 3)) %>%
  reactable()

We need to make sure and save our output file(s).

# Choose topTags destination
dir_output_limma <-
  path.expand("~/Desktop/Genomic_Data_Analysis/Analysis/limma/")
if (!dir.exists(dir_output_limma)) {
  dir.create(dir_output_limma, recursive = TRUE)
}

# for shairng with others, the topTags output is convenient.
top.table %>% tibble() %>%
  arrange(desc(adj.P.Val)) %>%
  mutate(adj.P.Val = round(adj.P.Val, 2)) %>%
  mutate(across(where(is.numeric), signif, 3)) %>%
  write_tsv(., file = paste0(dir_output_limma, "yeast_topTags_limma.tsv"))

# for subsequent analysis, let's save the res object as an R data object.
saveRDS(object = res, file = paste0(dir_output_limma, "yeast_res_limma.Rds"))

Be sure to knit this file into a pdf or html file once you’re finished.

System information for reproducibility:

pander::pander(sessionInfo())

R version 4.2.2 (2022-10-31)

Platform: aarch64-apple-darwin20 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats4, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: org.Sc.sgd.db(v.3.16.0), AnnotationDbi(v.1.60.2), IRanges(v.2.32.0), S4Vectors(v.0.36.2), Biobase(v.2.58.0), BiocGenerics(v.0.44.0), edgeR(v.3.40.2), limma(v.3.54.2), reactable(v.0.4.4), statmod(v.1.5.0), BiocManager(v.1.30.21.1), pander(v.0.6.5), knitr(v.1.43), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.2), purrr(v.1.0.1), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.2), tidyverse(v.2.0.0) and pacman(v.0.5.1)

loaded via a namespace (and not attached): httr(v.1.4.6), sass(v.0.4.7), vroom(v.1.6.3), bit64(v.4.0.5), jsonlite(v.1.8.7), bslib(v.0.5.0), highr(v.0.10), blob(v.1.2.4), GenomeInfoDbData(v.1.2.9), yaml(v.2.3.7), pillar(v.1.9.0), RSQLite(v.2.3.1), lattice(v.0.21-8), glue(v.1.6.2), digest(v.0.6.33), XVector(v.0.38.0), colorspace(v.2.1-0), htmltools(v.0.5.5), reactR(v.0.4.4), pkgconfig(v.2.0.3), zlibbioc(v.1.44.0), scales(v.1.2.1), tzdb(v.0.4.0), timechange(v.0.2.0), KEGGREST(v.1.38.0), generics(v.0.1.3), ellipsis(v.0.3.2), cachem(v.1.0.8), withr(v.2.5.0), cli(v.3.6.1), magrittr(v.2.0.3), crayon(v.1.5.2), memoise(v.2.0.1), evaluate(v.0.21), fansi(v.1.0.4), tools(v.4.2.2), hms(v.1.1.3), lifecycle(v.1.0.3), munsell(v.0.5.0), locfit(v.1.5-9.8), Biostrings(v.2.66.0), compiler(v.4.2.2), jquerylib(v.0.1.4), GenomeInfoDb(v.1.34.9), rlang(v.1.1.1), grid(v.4.2.2), RCurl(v.1.98-1.12), rstudioapi(v.0.15.0), htmlwidgets(v.1.6.2), crosstalk(v.1.2.0), bitops(v.1.0-7), rmarkdown(v.2.23), gtable(v.0.3.3), DBI(v.1.1.3), R6(v.2.5.1), fastmap(v.1.1.1), bit(v.4.0.5), utf8(v.1.2.3), stringi(v.1.7.12), parallel(v.4.2.2), Rcpp(v.1.0.11), vctrs(v.0.6.3), png(v.0.1-8), tidyselect(v.1.2.0) and xfun(v.0.39)